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ABSTRACT: We calculate the transverse momentum dependence in the production of two 
back-to-back hadrons in electron-positron annihilations at the medium/large energy scales 
of Bes-III and Belle experiments. We use the parameters of the transverse-momentum- 
dependent (TMD) fragmentation functions that were recently extracted from the semi- 
inclusive deep-inelastic-scattering multiplicities at low energy from Hermes . TMD evo¬ 
lution is applied according to different approaches and using different parameters for the 
nonperturbative part of the evolution kernel, thus exploring the sensitivity of our results to 
these different choices and to the flavor dependence of parton fragmentation functions. We 
discuss how experimental measurements could discriminate among the various scenarios. 
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1 Introduction 

Transverse momentum dependent (TMD) parton distribution functions (PDFs) and frag¬ 
mentation functions (FFs) depend on the longitudinal and transverse components of the 
momentum of partons with respect to the parent hadron momentum, as well as on their 
flavor and polarization state. The TMD PDFs and TMD FFs enlarge the amount of nonper¬ 
turbative information carried by ordinary integrated PDFs and FFs because they open the 
window on explorations of the multi-dimensional structure of hadrons in momentum space 
in terms of their QCD elementary constituents. For example, in the last years several data 
for single- and double-spin asymmetries in semi-inclusive deep-inelastic scattering (SIDIS) 
have been accumulated and can be interpreted as originating from the effect of specific 
combinations of (polarized) TMD PDFs and TMD FFs (for a review, see Refs. [1-4]). 

The TMD PDFs and TMD FFs can be defined only by a careful selection of physical 
observables that are sensitive to processes with two separate scales. In addition one needs 
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to study the appropriate factorization theorems for these observables. For example, the 
appropriate factorization theorem for SIDIS holds true if the hard photon virtuality is 
accompanied by transverse momenta of the order of nucleon mass [5, 6], which then are 
observed as a mismatch of collinear momenta. It is necessary that the definition of TMD 
functions includes all factorizable long-distance contributions to the physical cross section. 
These nonperturbative contributions, related to collinear gluon radiation, are summed into 
socalled gauge links that make the TMD functions color gauge invariant objects. Gauge 
links provide also the necessary phase to generate the above mentioned spin asymmetries [7— 
9]. Because initial-state and final-state gluon interactions are summed into different gauge 
links, the TMD functions may be process dependent, although parity and time-reversal 
invariance can simplify this non-universality to a simple proportionality factor [10, 11]. To 
account for scale dependence, the TMD functions obey evolution equations that generalize 
the standard Renormalization Group Evolution (RGE) to a multi-scale regime in hard 
processes. TMD evolution equations have been derived for unpolarized TMD PDFs and 
TMD FFs [12, 13], and for polarized ones only in a limited number of cases [14-16]. But 
despite these recent achievements, the phenomenological implementation of these effects 
is still under active debate [16-19]. From the experimental point of view, only few data 
sets are available with enough statistics that allows for a multidimensional analysis and 
a direct access to transverse momentum distributions [20, 21]; in other cases, the studies 
were limited in the multidimensional coverage and by the restricted variety of targets and 
final-state hadrons [22-26]. 

In a preceding paper [27], the dependence of the intrinsic transverse-momentum distri¬ 
bution of both unpolarized TMD PDFs and TMD FFs upon the flavor and the longitudinal 
momentum of the parton involved was discussed using the recently published data from 
the Hermes collaboration [20] on multiplicities for pions and kaons produced in SIDIS off 
proton and deuteron targets. Although the flavor-independent fit of the data was not statis¬ 
tically excluded, a clear indication was found that different quark flavors produce different 
transverse-momentum distributions of final hadrons, especially when comparing different 
species of final hadrons. This feature corresponds quite naturally to the well known strong 
flavor dependence of integrated PDFs [28-31], and to indications from some models [32-37] 
and lattice calculations of TMD objects [38]. The SIDIS process is useful because it gives 
simultaneous access to TMD PDFs and TMD FFs. But the factorized cross section always 
involves a convolution of transverse momenta of the initial and the fragmenting partons: 
anticorrelation hinders a separate investigation of the two intrinsic distributions. Moreover, 
the Hermes data were collected at such a limited range in the hard scale that the statistical 
analysis of Ref. [27] was reasonably performed even without involving modifications due to 
evolution effects. 

In this paper, we consider the semi-inclusive production of two back-to-back hadrons 
in electron-positron annihilations. In analogy with the SIDIS process, we define the multi¬ 
plicities in e + e~ annihilations as the differential number of back-to-back pairs of hadrons 
produced per corresponding single-hadron production. Then, we study their transverse 
momentum distribution at large values of the center-of-mass (cm) energy, starting from an 
input expression for TMD FFs taken from the analysis of Hermes SIDIS multiplicities at 
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low energy performed in Ref. [27]. In this framework, we can extract clean and uncon¬ 
taminated details on the transverse-momentum dependence of the unpolarized TMD FF, 
which is a fundamental ingredient of any spin asymmetry in SIDIS and, therefore, it affects 
the extraction also of polarized TMD distributions. Moreover, we can make realistic tests 
on the sensitivity to various implementations of TMD evolution available in the literature, 
since the hard scales involved in e + e - annihilations are much larger than the average values 
explored in SIDIS by Hermes , which is assumed as the starting reference scale. 

An important difference between PDFs and FFs is the role of the gauge links arising 
mostly from resummation of gluons with collinear polarizations. T-odd effects for PDFs 
enter through the operator definitions of the PDFs after inclusion of appropriate gauge links 
having also transverse pieces. For FFs T-odd effects are contained in the hadronic states 
and as a consequence there are less universality-breaking effects for FFs [6, 39-41]. 

In Sec. 2, we outline the theoretical tools needed to work out the cross sections for 
annihilations in two hadrons and define the e + e _ multiplicities. In Sec. 3, we introduce 
the QCD evolution of TMD FFs as the action of an evolution operator on input fragmenta¬ 
tion functions, we describe some procedures to separate perturbative from nonperturbative 
domains of transverse momenta, and we provide some prescriptions to parametrize the 
nonperturbative contributions to the evolution kernel and the resummation of soft gluon 
radiation. In Sec. 4, we introduce the flavor decomposition of fragmentation processes. In 
Sec. 5, we make predictions for the spectrum in transverse momentum of e + e~ multiplic¬ 
ities for production of two back-to-back hadrons, focusing on the sensitivity of results to 
the flavor of the fragmenting parton and to the different prescriptions for describing TMD 
evolution. Final comments and remarks are summarized in Sec. 6. 

2 Multiplicities for e + e _ annihilation into two hadrons 



Figure 1. Kinematics for the e + e annihilation leading to two back-to-back hadrons with momenta 
Pi and P 2 . 


We consider the process e + e —> hih^X depicted in Fig. 1. An electron e and a 
positron e + annihilate producing a vector boson with time-like momentum transfer q 2 = 
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Q 2 > 0. A quark and an antiquark are then emitted, each one fragmenting into a residual 
jet containing a leading hadron that for simplicity we will consider unpolarized: the hadron 
hi with momentum and mass Pi, Mi, and the hadron h- 2 with momentum and mass P 2 , M 2 . 
The two hadrons belong to two back-to-back jets, i. e. we have Pi • P 2 ~ Q 2 . In the following, 
we will limit Q 2 values to a range where the vector boson can be safely identified with a 
virtual photon. Using the standard notations for the light-cone components of a 4-vector, 
we define the following invariants 


2 Pi • q _ P\ _ -Pi ' P2 _ 2P2 • q _ P 2 + _ P-2 ■ Pi 
Q 2 q-P 2 32 “ Q 2 ~ ~q+ ~ q-Pi 


y = 


P 2 -t 
P- 2 -q ’ 


( 2 . 1 ) 


where £ is the electron momentum. The z\ is the fraction of parton momentum carried by 
the hadron hi, and similarly for z 2 referred to the hadron h 2 . Covariantly, we can define 
the normalized time-like and space-like directions 


t 11 = j- and F = 


= ^ pi j _ 

Pl-q 2 z 2 Q 2 


Correspondingly, we can dehne the projector into the space orthogonal to z and t: 


P£q" + q^PZ 


(P 1 ' = g^ - tH u + z»z v = g liu - ±2±-L* ^ + O I 

P 2 ■ q \Q 2 


M 2 


( 2 . 2 ) 


(2.3) 


The lepton momentum is then given by 


^ + (y~\) + QVyi 1 -v) L > ( 2 - 4 ) 

where £j_ = L/\£±\ and £j_ = g^£ v . 

The g^ projects onto the space orthogonal to q and P 2 . The projector onto the space 
orthogonal to Pi and P 2 , namely in the hadron cm frame where Pi and P 2 have no transverse 
components, is given by 


9t = 9^ 


= g? + 


p[ J p» + p$p" 

Pi-P 2 

P£ft + q£P 2 


„ , M 2 

+ ° [ w 


P 2 ■ q 


where the non-collinearity is defined as 


q£ = <? - ^ 

zi 


DM 
1 2 

— = 9t 9u 
z 2 


PL 

Zi 


+ 0 l-r 


M 2 

Q 2 


. /O / ILV-Plu 

+ °[ W )=-s1— + ° 


m 2 \ 

~W) 


(2.5) 


( 2 . 6 ) 


In the electron-positron cm frame of Fig. 1, we dehne the angle 9 = arccos(£ • z/\£\) 
where z = —P 2 . It is related to the invariant y ~ (1 + cos#)/2. In analogy to the Trento 
conventions [42], we dehne the azimuthal angle 


cos (f) 


P 2 X £ Pl_l X P 2 
P 2 x £\ \P\i_ x P 2 
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so that Pi = (0, |-Pl_l| coscf>, |-Pi±| sin (f>, 0) in this frame, and in any frame obtained from 
this one by a boost along z. In general, the covariant definition is cos 0 = — qx ■ £±/\qx\- 
The cross section for the e + e _ annihilation into back-to-back pairs of unpolarized 
hadrons can be written in a factorized formula at low transverse momenta [12, 16, 43, 44]: 

da hlfl2 67ra 2 


dzi dz 2 dqdy Q 2 


A(y)P(Q 2 ,n) 


x E e ' 


db T b T J 0 (q T b T ) z{ Df hl (zi, b T ; Cl, A*) zf {z 2 ,b T ; C2, A-0 + {q ^ q) 


*2/ 


+ Y{ql/Q 2 ) + 0{M 2 /Q 2 


( 2 . 8 ) 


where qx = \qt\ and A (y) = | — y + y 2 . The P is the hard annihilation part. The 
D\{z, b T ', (, n) is the TMD FF in impact parameter space for an unpolarized quark 
with flavor q fragmenting into an unpolarized hadron h and carrying light-cone momentum 
fraction 2 and transverse momentum conjugated to bx [45]. Both P and D q ~“ h are separated 
at the renormalization/factorization scale fj, and evolve with it through renormalization 
group equations. The D q ^ h depends also on the scale C (with ^ 1^2 = Q 4 ) aR d evolves with 
it via a process-independent soft factor. The term Y(q^/Q 2 ) ensures the matching with 
perturbative calculations at large transverse momenta. 

In this paper, we will consider a kinematics where q^ <C Q 2 and M 2 <C Q 2 ■ Hence, in 
Eq. (2.8) the Y(q^/Q 2 ) term and corrections from higher twists of order M 2 /Q 2 or higher 
will be neglected. Moreover, the soft gluon radiation is here resummed into the TMD 
FF at the Next-to-Leading-Log level (NLL). It implies that the hard annihilation part is 
consistently calculated at leading order (LO) in a s , namely P(Q 2 ,/j 1 ) ~ 1. Equation (2.8) 
then simplifies to 


da hlh2 


67rcr 


dz\ dz 2 dqj, dy Q 2 


A (y) 


E 


dbx bx Jo(qxbx) 


z\ D'f hl (zi, bx; Cl, d) 4 Dr h2 (z 2 , bx; C2, y) + (q ^ q ) 

(2.9) 


In Sec. 5, we present our results for the qx spectrum of hadron pair multiplicities 
in e + e _ annihilation. In strict analogy with the SIDIS definition [20], we construct the 
e + e~ multiplicities as the differential number of back-to-back pairs of hadrons produced 
per corresponding single-hadron production after the e + e~ annihilation. In terms of cross 
sections, we have 


M h ^(z 1 ,Z 2 ,q 2 T ,y) = 


da hlh2 


/ 


do h1 


( 2 . 10 ) 


dz\ dz 2 dqj, dy ‘ dz\ dy 

where da hlh2 is the differential cross section of Eq. (2.9). The da hl describes the production 
of a single hadron h\ from the e + e~ annihilation and it is obtained from the previous cross 
section by summing over all hadrons produced in one emisphere [43]: 

da hl 127ra 2 


dz\dy Q 2 


A (y ) 


<? 


e-Df- 


l (^i). 


( 2 . 11 ) 
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3 TMD evolution of fragmentation functions 


In the following, we describe in more detail the dependence of the fragmentation functions 
Df h of Eq. (2.9) upon the renormalization/factorization scale /r and the scale (. Different 
scenarios are possible according to the choice of the initial starting value for the factorization 
scale, and of the low-energy model describing the nonperturbative part of the evolution 
kernel. We first describe the structure of the input D q ^ h at the starting scale. 


3.1 Input fragmentation functions at the starting scale 

We consider the unpolarized TMD FF extracted by fitting the hadron multiplicities in 
SIDIS data at low energy from Hermes [20]. The assumed functional form displays a 
transverse-momentum dependent part which is described in impact parameter space by the 
following fixed-scale flavor-dependent Gaussian ansatz 1 : 


D^ h (z,b T ; Q 2 ) = df h (z; Q 2 ) exp 


a->ht 


i 

4 z 2 


[pir h (z)b 


(3.1) 


where ( Pj_) a h {z) with a = q, q. is the flavor- and ^-dependent Gaussian width at some 
starting scale Qq [27, 46, 47]. The choice of having separate Gaussian functions for different 
flavors is motivated by the significant differences displayed by the Hermes data between 
pion and kaon final-state hadrons [20]. The factorized collinear dependent part Q 2 ) 

is described by using the DSS parametrization of Ref. [48]. 

Following Refs. [49, 50], a possible energy dependence of the Gaussian distribution was 
taken into account introducing the logarithmic term 

e* P {-®^1.7|}, (3.2) 

with g 2 a free parameter. Choosing Qq = 1 GeV 2 , it was soon realized that the best-fit value 
for <72 w as compatible with zero. As a matter of fact, the Q 2 range spanned by Hermes 
is small and the obtained experimental data for multiplicities are not sensitive to evolution 
effects. For this reason, the fit was performed by using Eq. (3.1) at a scale fixed to the 
experimental average value, namely Q 2 = Q\ = 2.4 GeV 2 . With this choice, the possible 
energy dependence of Eq. (3.2) is automatically eliminated. 

In summary, the input to our studies on the evolution of Df^ h with the scales /r and 
£ is referred to the expression in Eq. (3.1) to be considered at the starting scale Q q = 2.4 
GeV 2 . However, depending on the choice of the initial value of the factorization scale this 
identification is not always straightforward, as will be explained in the following sections. 


3.2 The fib prescription 

As shown in Eq. (2.9), the TMD FFs generally depend on the factorization scale fi and 
on the scale £, that for convenience we name the rapidity scale. The TMD FFs satisfy 

1 The 1/z 2 factors appearing in Eq. (3.1) are due to br being conjugated to the partonic transverse 

momentum kr, whereas the TMD FFs in Ref. [27] are defined and normalized in momentum space with 
respect to the hadronic transverse momentum Kt = —zkr- 
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evolution equations with respect to both of them [12, 13]. The evolution with respect to 
is determined by standard RGE equations, whereas the evolution in ( is determined by a 
process-independent soft factor [12, 13]. 

The functional form of TMD FFs at small bx can be calculated in perturbative QCD. 
Conversely, the nonperturbative part at large bx must be constrained by fitting experi¬ 
mental data. At the medium/large energies of the Bes-III and Belle experiments, the 
perturbative tail of TMD FFs needs to be taken into account. Using the technique of 
Operator Product Expansion (OPE), it can be represented as a convolution of (perturba- 
tively calculable) Wilson coefficients C with the (nonperturbative) collinear fragmentation 
functions d± [12, 13]: 

D a - h (z,b T -, C,M) = lC®dr h ](z,b T -X^) + 0(bT Aqcd) • (3.3) 

'-V-' V - v/ -' 

small bx large bx 


The convolution is defined as 


[C®dr h ]{z,b T -, C,n)= [ — Cj-n (~,b T ;C,lj) (s\ n) . 

- J z $ \ S / 


(3.4) 


j=q,<hg 


The dependence of the coefficients upon both factorization and rapidity scales can be rep¬ 
resented in a factorized form: 


Cj^a(z,b T -, £,n) 



—K(b T ',g) 


Cj^a(Zibx]nl,n) , 


(3.5) 


where is defined as 




2e~ lE 

bx 


(3.6) 


and 7 e is the Euler constant. The K function in Eq. (3.5) 2 arises from the process- 
independent soft factor that is necessary to proof the factorization theorem leading to the 
definition of the TMD FFs; it drives the evolution of TMD FFs in the ( variable. The 
convolution in Eq. (3.4) is only valid for small bx, namely bx <C 1/Aqcd- Moreover, the 
expression of the C coefficients consists in a power series in a s In (/r 2 /V 2 ) (including also 
double logarithms of the same argument). The OPE is valid only when the logarithms 
do not diverge; this is accomplished, e.g., by choosing [i = ^ or qx, so that the series 
converges. Accordingly, if we choose /i = /ij we can write the TMD FF as 


D a - h (z,b T ;C,db) 



—K(b T ;ii b ) 



ds 

S 


c . 





+ 0(bx Aqcd) • 


d'C h (s- n b ) 

(3.7) 


The evolution of this fragmentation function from fi b to another value of /r (e.g., n = Q) 
is driven by RGE equations. Instead, the evolution from an initial rapidity scale Q to £ is 

2 Our K function corresponds to the D function in Ref. [13], and to the I\ function in Ref. [12] but for 
a factor —1/2. 
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controlled by the K function. The final expression of the TMD FF at the scales /x = Q and 

C is 


D a ~ h (z,b T ](,Q) = exp 


f^ =Q dp, 

/ —Iff 
),M, V 
—K(b T ;fi b ) 

»V £ 


C 
Ci 

ds 


—K(b T ;fi b ) 


Cj->a ( ibT] fJ^bi f^bj {s;fl b ) 


j=q,q,g 


+ 0(br Aqcd) , 


(3.8) 


where the anomalous dimension 7 ff reads 

IFF = ~ ^r cusp ln-^ + 7 V ^j , (3.9) 

and r cusp and 7 ^ are also power series in ot s in the MS scheme [16]. 

The above procedure is valid up to a maximum value of bx, that we name b max , beyond 
which we do not trust the perturbative calculation. Hence, it is convenient to reconsider 
the OPE by introducing the new variable bx that freezes at 6 max when bx becomes large: 

lim MM=Vax- (3.10) 

67 1 —1-00 


For bx iS frmax; the evolution in £ is controlled by the function K(bx ; /xg), where 




2e~ lE 

bx 


(3.11) 


The nonperturbative part at large bx is defined as what is left over [51]: 


9n P {bx) = -I<{bx;9b) + K ( b T',Vb) ■ 


(3.12) 


By adding the intrinsic transverse distribution at the starting scale (see Eq. (3.1)), 
Eq. (3.8) becomes 


D a ~ h {z, b T ; C, Q) = ex P | J ^ 


—K(b T ;ii b )-gnp(b T ) 


Vi 


£ 


- K(b T ;p b )-gn P (bT) 


Sr;/xf.Ms) d V h ( s iH) 


3=9,9,9 


x 1 

~ \Q 2 o 


-9np(b T ) 


(3.13) 
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If we insert Q = /if and £ = fi 2 = Q 2 , the above equation reduces to 


D a " h (z,b T ;Q 2 ,Q) = ex: pi f ^ 7 FF 

I V 


g_ 

& 


2 \ - K {i>T;ni)-gnp{b T ) 


E 

j=q,<h. 9 ' 


y Cj-v (^, S T ; Mg,Ms) ^(s;^) 


1 t 2 //if \ 9np(bT) 

x e 4*2 hi 0 \ 


z* 


Ql 


= R(b T ; Q , Q, 4 fi h ) D a ^(z, b T -lit, H ) . 


(3.14) 


Hence, the net effect of evolution can be represented as the action of an evolution operator 
R on the input TMD FF evaluated at the scale /ig, which is running with b F - This peculiar 
feature grants that there is a smooth matching between the perturbative domain at small 
bT and the nonperturbative domain at large b F - It is interesting to remark that from 
Eqs. (3.13) and (3.14) we deduce that modelling the nonperturbative part affects the whole 
b F spectrum, not only the large b F region. 

In this paper, we resum the soft gluon radiation up to NLL contributions in In (////!&), 
which corresponds to include terms linear in ot s in the perturbative expansion of K and 7 ' , 
and quadratic in the expansion of r cugp [16]: 


, C F a 2 
A(6 t ;^) = — a s In— , 
2 tt nf 


1 V = 


3 Cf 


OLr 


r - Cf fi 1 as 
1 cusp - Q7 S 1 T ^ 


67 

~9~ 




(3.15) 


where C A = N c , C F = (. N, 2 - l)/2 N c , are the usual Casimir operators for the gluon and 
fermion representations of the color group SU(1V C ) with N c colors, and T F = rif/2 with nj 
the number of active quark flavors. Consistently, the coefficients C are computed at LO in 
ac s , namely they reduce to S functions such that Eq. (3.14) simplifies to 


D a - h (z,b T ]Q) = exp 


rQ dp, 

—1 


FF 



—^NLL (Pt ; A7 ) — g np (&T ) 


X dT ( 2 ; Mg) 72 eX P { - 


(pir h r) 

4 z 2 



-9np(b T ) 


(3.16) 


The definition of g n p(bx) in Eq. (3.12) obviously implies that this function depends 
on 6 max , i.e. on the value of the impact parameter that sets the separation between the 
perturbative and nonperturbative regimes. Indeed, by perturbatively expanding K(b F ; fJ-b) 
at lowest order we have [51] 

S„ P (M«^i)^ ln ( 1+ . (3.17) 

F V ®max J 
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For bx <C 6 max , this expression recovers the quadratic parametrization \g2^>x adopted in the 
fits of Refs. [50] and [52], and it suggests that the parameter g 2 is not free but anticorrelated 
to 6 max , and proportional to 6^ ax through a perturbatively calculable coefficient. The g np 
function accounts for the radiation of soft gluons emitted from a parton. A small (large) 
value of & max implies that the QCD perturbative description is valid up to relatively small 
(large) bx values. Consequently, the amount of soft gluons emission is larger (smaller) and 
we expect a large (small) value for 52 • More generally, this anticorrelation is motivated by 
the fact that both the exact function K(bx', gb) and the TMD FF itself must not depend on 
the arbitrary choice of 6 max . So, 6 max should not be regarded as a free parameter to be fitted 
to data, but it should be considered as an arbitrary scale that separates perturbative from 
nonperturbative regimes: changing 6 max implies a rearrangement of all terms in Eq. (3.13) 
such that the TMD FF does not change [51]. 

For the purpose of this work, we will consider anticorrelated pairs of values for {6 max , 52 }) 
inspired to the values adopted in Refs. [50] and [52], We will also explore different expres¬ 
sions for each one of the bx and g np functions. For bx, our first choice is the socalled “6-star” 
prescription [12, 50] 

bx 


b T = b* T = 


1 + 


by 


(3.18) 


The second choice is based on the exponential function 


Lt — bji — &max ^ 1 exp 


b 4 

u maxJ 


(3.19) 


that is steeper and it approaches the asymptotic constant 6 max more quickly. For g np , we 
choose a linear function of b\ similarly to Refs. [49, 50, 52] (see also Eq. (3.2)): 


_lm/i \ _ 92 ,2 
Snpi^T) ~ b T . 


(3.20) 


The second choice is suggested by Eq. (3.17): 


5n P g (M = 92 In ( 1 + h ) . 


(3.21) 


This expression was considered also in Ref. [53], and it reduces to Eq. (3.20) for small bx- 
In principle, we have four different combinations of prescriptions: {&?),, (/])“}, {b^-,g l np}, 
{ 6 ^, g^p}, and However, after some preliminary exploration we realized that 

some of them were producing redundant results. Therefore, they have been neglected. In 
summary, the transverse-momentum spectrum of the multiplicities in Eq. (2.10) will be 
analyzed by varying the anticorrelated pair of parameters {b max ,g 2 }, an d by considering 
only the two combinations { 6 ^, 5 ]]°} and {fej-,, 

Finally, we remark that if we choose Q = g~ b in Eq. (3.16), i.e. if we switch off evolution 
effects, we should recover the Gaussian model expression of Eq. (3.1) for the TMD FF at 
the initial scale Qp. Formally, this is not the case because in the second line the collinear 
d\ is evaluated at and the term {g^/Q‘p)~ 9 ' [iV ^ bT ^ survives. However, the Gaussian model 
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of Eq. (3.1) is deduced by fitting the Hermes SIDIS data, whose kinematics overlaps the 
domain of very large bx 6 max , namely where bx ~ b m ax - If we use the prescription bx = b^ 
of Eq. (3.18), it is easy to check that for 6 max = 0.7 GeV -1 we have nf ~ Qo = 2.4 GeV 2 . 
Hence, the D a ^ h (z, bx', H%) of Eq. (3.16) at Q = /i^ actually behaves like the D a ^ h (z, bx', Qo) 
of Eq. (3.1) at the scale Qo and at very large bx values, or equivalently for very small parton 
transverse momenta. 


3.3 The fixed-scale prescription 

In Eq. (3.14), we have expressed the evolved TMD FF at a scale Q as the result of an 
evolution operator R acting on the same TMD FF evaluated at the scale Hi running with 
bx- Alternatively, we can fix the initial scale at the value Q? = Qq = 2.4 GeV 2 for the 
whole bx distribution: 


D°- h {z, b T ; Q) = R(b T ; Q, Qi) D a - h {z, b T ; Qi) - (3.22) 

With this choice, it is not possible to apply the OPE for calculating a perturbative tail to 
which the TMD FF should match at low bx, as it was done in Eq. (3.3): we need a model 
input over the whole bx spectrum. In our case, it is now very easy to identify the input 
TMD FF at the starting scale Qi with the Gaussian parametrization of Eq. (3.1) at Qo- 
Then, for p 2 = Q = Q 2 = Qq = 2.4 GeV 2 the TMD FF evolved at NLL up to a final scale 
H 2 = C = Q 2 becomes 


D a ~* h (z, b T ; Q) = exp { / ^7 ff 




'Qi V 


I NLL 


2 \ —-KnLI X b T',Qi) 


xdr h iz-Qi) 


/ p2 \ a ~*h 


exp 


(*) 


4^ 2 


b% 


(3.23) 


The contribution from the g np term in the input distribution does not appear because of 
the choice of the starting scale Q = Qj = Q q. 

The choice Hi = Qi of identifying the starting factorization scale with a fixed scale for 
the whole bx spectrum has important consequences also on the function K. From Eq. (3.15), 
we can expand K in powers of In (n/^b)'- if Hi 7 ^ db, the series may not converge. One 
possible workaround is to apply the resummation technique to the K function itself [13]. 
Here, we will discuss two different prescriptions: computing K from Ref. [12] at a fixed order 
in or dressing K by resumming large logarithms of the kind In (n/d'h) [13]. In the first 
case, K is expanded in powers of a s ; in the second case, the expansion is in a s In (n/^b)- 
If Hi = Hb, th e t wo expansions are the same. 

We will refer to the first choice as the "fixed-scale" prescription. Contrary to the 
prescription described in the previous section, there is no need to define an arbitrary scale 
iimax to separate perturbative from nonperturbative regimes. However, the function I\ is 
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evolved from m to Qi through its anomalous dimension: 


rQi du 

K(br ; Qi) = K(br', A*g) + / — r cusp + g np (bx ) 

•'ms M 

« rQi du 

NLL / T cusp + g np (b T ) , (3.24) 

where g n p(bx) can get either the expression in Eq. (3.20) or in Eq. (3.21). The pertur¬ 
bative contributions are calculated at NLL as in Eq. (3.15), according to which we have 
-Knll(&t; db) = 0. 

The second choice is connected to the results of Ref. [13], because we resum all large 
logarithms of the kind In (d/d'b) in the perturbative part as 

K{br\ Qi) = D R {b T ; Qi) 9(bx, c ~ bx) + gnp(br) 0(bx - b T , c ) , (3.25) 

where D R is the resummed contribution computed in Ref. [13], and bx, c is the convergence 
radius of the perturbative expression. Apart from the resummation of logarithms, the 
main difference with Eq. (3.24) is the presence of the 9 functions: no bx prescription is 
used to connect the perturbative and nonperturbative domains. And the nonperturbative 
contribution acts differently: while g np in Eq. (3.24) applies to the whole bx spectrum, in 
Eq.(3.25) it does only for bx > bx , c • Hence, we use the notation g np to account for this 
difference. For example, the K function must be at least continuous at bx = bx, c ■ We can 
match this constraint by defining the nonperturbative contribution at bx > bx, c as 


gn P (bx) = D R (bx, c ) 


1 + 9n P (bx ~ b T ,c) > 


(3.26) 


where g np can be again either the g^ prescription of Eq. (3.20) or the gn P prescription of 
Eq. (3.21). 

For Qi <C Q, the perturbative component D R in Eq. (3.25) diverges for bx < bx, c - 
Hence, its contribution to the evolution of the fragmentation function becomes negligible, 
being of the kind ( Q/Qi)~ D ■ Since K is a smooth function in bx, also the contribution of 
the nonperturbative part g np for bx > bx,c becomes numerically negligible [13]. However, 
this result cannot be generalized to any value of Q. Since we will make explorative calcu¬ 
lations also at the Bes-III scale Q = \/l4.6 GeV which cannot be considered to be much 
larger than the initial scale Q o = \/2.4 GeV of our input TMD FF, we will consider only 
the "fixed-scale" prescription of Eq. (3.24). 

3.4 Summary of evolution kernels 

In summary, we consider two possible ways of evolving the TMD FF, according to the 
choice of the initial factorization scale p,,;. It is understood that all formulae are computed 
at the NLL level of accuracy, according to Eq. (3.15). 
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A The "Hb" prescription: then g 2 = gf = Q, £ = p 2 = Q 2 , and we have 


£> a ^0, b T ; Q) = exp <( / )> I ^ 


dfi 


Q 


2 \ - K (.b T ;v- b )-gnp(b T ) 


'H 




x dr h ( 




exp 


(pir h (z) 

\z 2 


fry 


A 

4 


2 \ —9np(b T ) 


(3.27) 


where 7 ff and /i are described by Eqs. (3.9) and (3.15), /rg is given by Eq. (3.11) 
with Eqs. (3.18) and (3.19), and g np is described in Eqs. (3.20) and (3.21). 

B The "fixed-scale" prescription: then /r 2 = Qf = Q, f = g 2 = Q 2 , and we have 


D a ~ (. z , b T ; Q) = exp <J j —7 ff j y Q? 


, ' <3 dg 


Q 2 ^-K(p T -,/i i )-f?i fr cusp -g np {b T ) 


x di~" h (z; Qi) A- exp <( - 


1 


^(s) , 

4z 2 


(3.28) 


where 7 fFi Mg; Snp are dehned in the same equations as above, while it is given in 
Eq. (3.24). 


4 Flavor dependence of fragmentation functions 

The flavor sum in Eq. (2.9) can be made explicit and further simplified using the symmetry 
upon charge-conjugation transformations: 

Dr h (z, b T ; Q 2 ) = Df*(z, & T ; Q 2 ) . (4.1) 

At the starting scale Q q, we distinguish the favored fragmentation where the fragmenting 
parton is in the valence content of the final hadron h. All the other channels are classified 
as unfavored fragmentation and are characterized by the fact that the detected hadron is 
produced by exciting more than one qq pair from the vacuum. If the final hadron is a kaon, 
we further distinguish a favored fragmentation initiated by an up quark/antiquark from the 
one initiated by a strange quark/antiquark. We limit the sum to three flavors u, d, s, and 
the corresponding antiquark partners. 


4.1 Favored and unfavored fragmentation to different hadron species 

For the final hadron pair being (hi, h 2 ) = (vr + , 7r _ ), the flavor sum in Eq. (2.9) becomes 


E e * D T + Dr +(q^q) = Dir- + D 


y unf i 


(4.2) 


where 


Dir~(z 1 ,z 2 ,b T ;Q 2 o) = lDr + r,b T -,Q 2 o)Dr-(z 2 ,b T -,Q 2 o) 

+ \ D r + (*!, 6 t ; Ql) D d r~ (*2,6 t ; Ql) 


(4.3) 
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and 


D^r(zi,z 2 ,b T ;Q 2 0 ) = ^Df^ + (z 1 ,b T -,Q 2 o)D^-(z 2 ,b T ;Q 2 0 ) 

+ (zi,b T ; Ql) Dt**- (z 2 , b T - Ql) 

+ i (dp + (*!, 6 t ; gg) or'(* 2 , br, Ql) 

+Dr + (zi, b T ; g§) (* 2 , 6t; g§)) • (4.4) 


Using the charge-conjugation symmetry of Eq. (4.1), it is simple to prove that the result 
for (hi, / 12 ) = (tt - , 7r + ) is identical to the above one in Eq. (4.2). 

If the final pions have the same charge, (hi, / 12 ) = (vr + , vr + ), we have 


J 2 el D r +D r + +( q ^q) 


4 4 

-Dr* 

9 1 


(zi,br] Ql) D^ ni (z2, b T ; Ql) 


+ ^ D^(zi, b T ', Ql) Df^ + (z 2 ,b T - Ql) 


+ (!<-»• 2) 


+ g D lt f(^i. b T ; Qo) -°unf (- 2 , br ; g§), (4.5) 


where 


£>?&(*, &r;Qo) = ^ + («,6 t;Q§ ) = J Df- r+ (z,6 T ;Q2 ) 

= Z^ + Mt;Q^) = 44r + (z,6 r ;Qg) . (4.6) 

Again, because of charge-conjugation symmetry we get the same result for (hi, / 12 ) = 
(vr", TT~). 

If the final hadron pair is (hi, I 12 ) = ( K + , K~), the flavor sum becomes 

E °r K+ d V k ~ + q) = d£ k ~ + d£«- , (4.7) 

I 

where 

D£ K -(zi,Z2,b T -Ql) = ^Dr K+ (zi,b T -,Ql)Dr K ~(z 2 ,b T -,Ql) 

+ ^Dr K+ (zi,b T ;Ql)Dr K ~(z 2 ,b T -,Ql) , (4.8) 

and 

D^ K -(zi,z 2 ,b T -Ql) = ^Dr K+ (zi,b T ;Ql)Dr K ~(z2,b T -,Ql) 

+ ^Dr K+ (zi, b T ;Ql) D\- k ~ (Z 2 , b T ; Ql) 

+ l (Dt K+ (zi,b T - Ql) Dt K ~ (z 2 , br, Ql) 

+D d r K+ ( Zi , b T -Ql) Dt K ~ (Z 2 , br Ql)) • (4.9) 
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Charge-conjugation symmetry grants the same result for (hi, h 2 ) = (K , K + ). 
If (h!, h 2 ) = (K+, K+): 


D T K+ d V k+ + (q^q) = 


Dr K+ (zi,b T -Ql)D^(z 2 ,b T -Ql) 


2\ r > K + , 


9 


+ ^ DS (zi,b T - Ql) Dt K+ (z 2 , b T - Ql) 


+ (!<-)• 2) 


where 


+ g ^unf ( z l^T\Q q) £>unf ( z 2, ^Ti Qq) , (4-10) 


(z, b T ; Ql) = D r f K+ (z, b T - Ql) = Df K+ (z, b T ; Q§) 

= Dt K+ (z,b T ;Q 2 0 ) = Dt K+ (z,b T -Ql) . (4.11) 


As before, we get the same result for (hi, h 2 ) = ( K , K ). 

The last combination is (hi, h 2 ) = (tt + , K ~): 

£ e ' D i^ + D V K ~ + ?) = (*i, *2, 6 t; Q§) 

q 

+ ^ (z?f^ + ( ai , 6r; Qo) C(* 2 , Q§) + ^ n + f (z,, b T ■ Ql) DT K ~ (z 2 , b T ; Q§)) 

+ 2 -^unf^l) t>T', Qo) ( z 2, &T1 Qo) > (4-12) 

where 

^-(zi,z 2 ,6 T ;Qg) = l^- r+ (^ 1 ,6 T ;Q2) jD ^- (22)6T ;Q2 ) (4.1-3) 

and charge-conjugation symmetry applied to Eq. (4.11) gives 

D^(z,b T -Ql) = Dr K -(z,b T ;Ql) = Dr K -(z,b T ;Ql) 

= Dt K '(z,b T ;Ql) = Dt K ~(z,br-,Ql) , (4-14) 

and grants that the same result in Eq. (4.12) holds also for (hi, h 2 ) = (tt~, K + ). 


4.2 Flavor dependent Gaussian ansatz 

The starting input to our analysis are the TMD FFs extracted by fitting the hadron multi¬ 
plicities in SIDIS data from Hermes at Ql = 2.4 GeV 2 [27]. The assumed functional form 
displays a transverse-momentum dependent part which is described in impact parameter 
space by the following flavor-dependent Gaussian ansatz: 

= dr h (z-Ql)G h q (z,b 2 T ). (4.15) 

The cross section of Eq. (2.9) (and, in turn, the multiplicity in Eq. (2.10)) is then a sum of 
Gaussians, and thus no longer a simple Gaussian. The width of the Gaussian depends also 


D c r h (z, b T \ Ql) = d1~* h (z; Ql) exp 
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on the fractional momentum z, as done in several model calculations or phenomenological 
extractions [36, 54-58]. The chosen functional form is [27] 


(pir h ( Z )=(pi) 


2 \ q~"h (z? + (5) (1 - Zp 


(z@ + (5) (1 — zp 


(4.16) 


where j3, 5 , 7 , are fitting parameters and (P±) q h = ( P± ) q h (z), with z = 0.5. 

Isospin and charge-conjugation symmetries suggest four different Gaussian shapes [27]: 


(Piy^ = = (Pi)™~ = (Pl) d ^~ = (Pl) iav 

*K+ 


^ p 2 \ ~ T ' _ / 732 \ u~>K _ 


=(m 


(piy^ K+ =(pi)^ K ~ = 

all others _ 


_ (pjyK 

- (piy K ■ 


(4.17) 

(4.18) 

(4.19) 

(4.20) 


Correspondingly, we have four different Gaussian functions in Eq. (4.15): 


G TT' /^T7T ' S~17T S~ 11 T ( r2 \ 

u — — U u — = 'JfavfZj 0 T ) . 

G« + =G«~ =G uK (z,b 2 T ) 


Gf + = Gf = G sK (z,b^) 


iK~ _ 


G tt _ s-iK _ S~llT~ r _ _ s~nr _ r" 1 " _ n&' 

0 , — ^ rl — Lt c — Lt c — 


= Gj = Gf * = Gf = Gf = G un f(z, bp 


(4.21) 

(4.22) 

(4.23) 

(4.24) 


Each one of these four functions depends on the same /3, 6 , 7 , fitting parameters of Eq. (4.16), 
such that all the Gfz,bf in Eq. (4.15) are described by seven parameters. 

For the collinear functions df h (z\ Qq), we adopt the same assumptions of Ref. [48]: 


- isospin symmetry of the sea quarks 

- for h = 7 r + , a direct proportionality between the {d + d) and (u + u) combinations, 
he. (d + d) = N(u + u). 


Therefore, we have three independent favored fragmentations, 

d'Tfz-Ql), dT K \z-Ql ), df K+ (z-Ql), (4.25) 

and two independent unfavored fragmentations: 

<Cf* = if* = if* = if* = <Cf(z; Ql) , (4.26) 

<T* + = if K * = if** = if** - Ol Ql) ■ (4.27) 

The remaining favored channel d —> 7r + is then given by 

df* + = N d^ + + (N — 1) df^ + 

= N d?T + + {N- 1) df uni (z; Ql) . (4.28) 
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The h = n~ and h = K~ channels can be deduced from the above ones using the charge- 
conjugation symmetry of Eq. (4.1). 

By inserting the Gaussian ansatz with the above assumptions in the expressions of 
Sec. 4.1, we get 


Dl v ^(z 1 ,z 2 ,b T ;Q 2 0 ) = 


N +4 2 \ JU^7T+/„ . 


9 


-dr {zi-Q^dr (z 2 -Q 2 0 ) 


+ N{N q l) dr\zr,Ql)df u Az2-, Ql) + (l o 2 ) 

+ “ ~ —d* im f(-gi; Ql) d\ nn{ (z 2 ] Qo) Gf av (zi, bji) Gf av (z 2l bx) , 


9 


7 


D lli ( Z \, Z2, b T ; Ql) = ^ dJlAzi; Ql) d\ urii (z 2 \ Ql) G unf (zi, b 2 T ) G unf (z 2 , b\) , 
D u r\zM-,Ql) = dr + (z-,Ql)G fav (z,b 2 T ) , 

Dr + (z,b T -Ql)= [Ndr^z-QD + iN-^df^iz- Ql)] G iav (z,b 2 T ) , 
AT nf ( z , br', Ql) = di nn f(z; Ql) G un f (z, b , 

d (C K (zi, z 2 , b T ; Ql) = ^dr + (zi ; Qo) d 1 f K+ (z 2 - 1 Qo) GW(*i, &t) G uK (z 2 , 6^) 
+ ^df x+ (z i; Q 2 ) dr^ + (z 2 ; Qo) G a *(*i, 4) G sK (z 2 , b 2 T ) , 

A ( 21 , 22 , &t; Qo) = ^ rff u + nf(2i; Qo) di^ nf {z 2 -, Qo) G unf (zi, b 2 T ) G unf {z 2 , b |) , 

^r K+ (2, b T ; Qg) = dr A ' + (2; Qg) G«a-(2, #) , 

T»f K+ (z,6 T ;Q§) = £>r*>,&r;Qo) = dr* + (*;Qo)G s *Mr) , 

D un{ (z,b T -,Ql) = dfunf(2; Qo) G un f(z, , 


tf£*"( 2 i,*»,&r;Q 8 ) = ^dr^z 1 -Ql)dr^r,Ql)G { rzi,b 2 T )G uK (z 2 ,b 2 T ) . 


2 \ ru->K + / 


(4.29) 

(4.30) 

(4.31) 

(4.32) 

(4.33) 


(4.34) 

(4.35) 

(4.36) 

(4.37) 

(4.38) 

(4.39) 


5 Predictions for TMD multiplicities 

In this section, we present our results as normalized multiplicities 

M h ^( Zl , z 2 , q 2 T , y)/M h ^(z u z 2 , 0, y) (5.1) 

for the hadron pair (h\,h 2 ), where M hlfl2 (zi, z 2 ,q^,y) is defined in Eq. (2.10). In such 
way, we are able to directly compare the genuine trend in q^ for each different case. If 
not explicitly specified, we choose y = 0.2. For selected values of {z\. z 2 }, the results are 
displayed as a function of P 2 j_ = zfq Hence, the useful range in P 2 ± depends on z\ in 
order to fulfill the condition qj, <C Q 2 . The range obviously depends also on the choice of the 
hard scale; we consider Q 2 = 100 GeV 2 , as in the Belle experiment, and Q 2 = 14.6 GeV 2 , 
as in the Bes- 111 one. For each specific case, the results are displayed as uncertainty bands: 
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they represent the 68% of the envelope of 200 different values for the intrinsic parameters 
in Eqs. (4.16)-(4.20) for the D±(z, hr; Qq) at the starting scale Qg, obtained by rejecting 
the largest and lowest 16% of them. The 200 values are obtained by fitting 200 replicas of 
S1D1S multiplicities measured by the Hermes collaboration [20]. If the 200 values for each 
parameter were distributed as a Gaussian, the 68% band would correspond to the usual ler 
confidence interval (for more details, see Ref. [27]). 

The results are organized as follows. In Sec. 5.1, we show the sensitivity of the nor¬ 
malized multiplicity to different values of the evolution parameters {b max , < 72 } described in 
Sec. 3.2 for a final hadron pair (hi/ 12 ) = (7r + 7r - ). In Sec. 5.2, we compare normalized 
multiplicities for the two different evolution schemes described in Secs. 3.2 and 3.3. In 
Sec. 5.3, we discuss the capability of discriminating among the various prescriptions illus¬ 
trated in Sec. 3.2 for the nonperturbative evolution effects. In Sec. 5.4, we concentrate 
on the sensitivity of the normalized multiplicities upon varying the fractional energy 2 of 
final hadrons. In Sec. 5.5, we show how the results get modified when lowering Q 2 from 
the Belle scale to the Bes-III scale. Finally, in Sec. 5.6 we discuss the sensitivity of 
ratios of normalized multiplicities for different final states to the flavor structure of the 
intrinsic transverse-momentum-dependent part of the input TMD FF at the starting scale 
of evolution. 

5.1 Sensitivity to nonperturbative evolution parameters 

As already remarked in Sec. 3.2, for a specific evolution scheme the nonperturbative part 
of the TMD evolution depends on the choice of a prescription for describing the transition 
from perturbative to nonperturbative regimes, which in turn depends on the two parameters 
6 max and g 2 . In this section, we explore the sensitivity of our predictions to different values 
of the pair {6 max , < 72 }- We adopt as limiting cases the choices {6 max = 1-5, <72 = 0.18} and 
{ 6 max = 0.5, <72 = 0.68}, that were deduced in Refs. [52] and [50], respectively, by fitting 
the transverse-momentum distribution of lepton pairs produced in Drell-Yan processes. If 
not explicitly specified, the first choice is described by uncertainty bands with dot-dashed 
borders while the second choice is linked to bands with solid borders. As explained in 
Sec. 3.2, the two parameters are anticorrelated. In the following, we show results also for 
the interpolating choice { 6 max = 1, <72 = 0.43}. The corresponding results are displayed as 
uncertainty bands with dashed borders. 

In Fig. 2, the normalized multiplicity 

M 7r+7r_ (zi = 0.5, z 2 = 0.5, q%,y = 0.2) /M 7r+7r "(*i = 0.5, z 2 = 0.5,0 ,y = 0.2) (5.2) 

is shown as a function of P 2 ± = z 2 q = (0.5) 2 g^ at the Belle scale Q 2 = 100 GeV 2 for the 
"gb scale" evolution scheme and with the {h}, < 7 }'}} prescription for the transition to the 
nonperturbative regime, as explained in Sec. 3.2. The explored range in P 2 j_ is such that for 
z\ = 0.5 the maximum qj, satisfies the condition -C Q 2 ■ The three uncertainty bands, 
corresponding to the three different choices { 6 max = 1.5, <72 = 0.18} (dot-dashed borders), 
{ 6 max = 1, <?2 = 0.43} (dashed borders), and { 6 max = 0.5, (72 = 0.68} (solid borders), are 
well separated. The squared box with error bar indicates a hypothetical experimental error 
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Figure 2. The normalized multiplicity at Zi = Z 2 = 0.5 as a function of p l 2 x = zWt = (0.5 )Wt 
at the Belle scale Q 2 = 100 GeV 2 for the "/if, scale" evolution scheme and with the {6^, g} 1 "} 
prescription for the transition to the nonperturbative regime (see text). The uncertainty bands 
correspond to various choices of the nonperturbative parameters of evolution: {6 max = 1.5, 52 = 
0.18} for the band with dot-dashed borders, {6 max = 1, <72 = 0.43} for the one with dashed borders, 
{6 max = 0.5, g 2 = 0.68} for the one with solid borders. The latter is accompanied by a light-gray 
band with dot-dashed borders, that represents the result with the same parameters but with the 
choice /*;,/2 for the arbitrary matching scale, and by an overlapping light-gray band with dashed 
borders for the choice 2/if,. An experimental error of 7% is also indicated. 


of 7%. We fix it by propagating to the normalized multiplicity the typical experimental 
error of 3% for single-hadron production data in e + e~ annihilations at Q 2 = 100 GeV 2 and 
2 = 0.5, from which the collinear d\(z] Q 2 ) are extracted [48]. This experimental error of 
7% seems small enough to discriminate among predictions produced with different choices 
of {^maxj 92 }■ 

Two additional light-gray bands are shown, which are partially overlapped (dot-dashed 
borders) or completely overlapped (dashed borders) to the band with solid borders cor¬ 
responding to the choice {6 max = 0.5, 52 = 0.68}. These bands reproduce the outcome 
of calculations performed in the same conditions but for different (arbitrary) choices of 
the scale /x/,. If the band with solid borders corresponds to calculations with the choice 
of Eq. (3.6) for /if,, then the light-gray band with dot-dashed borders corresponds to the 
choice /Xf,/2, and the one with dashed borders to 2/x/. The almost complete overlap of these 
results shows that for the selected observable, the normalized multiplicity, the theoretical 
uncertainty in determining the matching scale /if, (that describes the transition from per¬ 
turbative to nonperturbative regimes) is negligible with respect to the sensitivity to the 
parameters describing the nonperturbative effects in the evolution. 

5.2 Sensitivity to evolution schemes 

In this section, we explore the sensitivity of our normalized multiplicity to the choice of 
the evolution scheme. In Sec. 3, we described two different schemes, the "/x/ scale" and the 
"fixed scale". They differ mainly in the fact that in the latter the whole distribution in 
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impact parameter space bx of the TMD FF D\ at beginning of evolution is computed at a 
fixed scale Qq, namely there is no impact parameter that describes the transition from low 
(perturbative) bx to high (nonperturbative) bx- Actually, one would expect that for small 
values of g 2 and corresponding not too large values of 6 max (he., where the perturbative 
description of the evolution of the bx distribution is still applicable and gives the predom¬ 
inant contribution) the predictions from the different schemes should tend to a common 
result, determined mainly by a fully perturbative calculation. However, the complexity 
of the evolution kernels, described in Secs. 3.2 and 3.3, indicates that this is too a naive 
expectation. 


M(Z!,z 2 ,Ph T ) / M(z lt z 2 ,0) 



Figure 3. The normalized multiplicity at Zi = Z 2 = 0.5 as a function of P'f_ L = zfq^ = (0.5 ) 2 q^ 
in the same conditions and with the same notation as in Fig. 2, but for the "fixed scale" evolution 
scheme. The additional light-gray bands with dot-dashed and solid borders are the result related 
to the "/i;, scale" evolution scheme for { 6 max = 1.5, gi = 0.18} and { 6 max = 0.5, gi = 0.68}, 
respectively. 

In fact, in Fig. 3 the normalized multiplicity of Eq. (5.2) is shown as a function of 
Pix = z \Qx = ( 0 . 5 ) 2 q^ at the Belle scale Q 2 = 100 GeV 2 with the { 6 ^, g^p} prescription. 
There are two groups of uncertainty bands. The former one displays the results for the 
"fixed scale" evolution scheme in the standard notation, i.e. for {6 max = 1-5, g 2 = 0.18} 
(dot-dashed borders), { 6 max = 1, <72 = 0.43} (dashed borders), and { 6 max = 0.5, g 2 = 0.68} 
(solid borders). Then, two additional light-gray bands are shown that correspond to the 
results with the "scale" evolution scheme for { 6 max = 1.5, g 2 = 0.18} (dot-dashed 
borders) and ( 6 max = 0.5, g 2 = 0.68} (solid borders). 

It is evident that for the maximum (minimum) 6 max ($ 2 ) the band with dot-dashed 
borders in the "fixed scale" scheme is not similar to the light-gray band with dot-dashed 
borders in the "scale" scheme. Actually, all the results in the "fixed scale" scheme 
show a much larger distribution in -P 2 ^, somewhat pointing to stronger evolution effects of 
perturbative origin that seem to be absent in the "^ scale" scheme (where the scale choice 
minimizes the effect of large logarithms in the perturbative coefficients). It is important to 
notice that there is a significant overlap between the band with dot-dashed borders in the 
"fixed scale" scheme and the light-gray band with solid borders in the "scale" scheme. 
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Apparently, the normalized multiplicity seems not to be enough sensitive to discriminate 
among different evolution schemes, since two different choices of them can produce similar 
results with different evolution parameters { 6 max , 92 }- However, this result is observed at 
a specific value of fractional energies of the final hadrons, namely Z\ = Z2 = 0.5. 


M(zi,z 2 ,Ph T > / M(z 1 ,z 2 ,0) 



Figure 4. The same as in the previous figure, but at z\ = 0.3. 

In Fig. 4, we show the P 2 ± distribution of normalized multiplicities calculated in the 
same conditions, notation and conventions as in the previous figure, but at z\ = 0.3 and 
Z 2 = 0.5. The band with dot-dashed borders in the "fixed scale" scheme can now be 
easily separated from the light-gray band with solid borders in the "fib scale" scheme if the 
indicated hypothetical experimental error is around 7%. Therefore, only when combining 
the study of both the z and P 2 ± dependencies in the normalized multiplicity we may be 
able to discriminate among different TMD evolution schemes. 

5.3 Sensitivity to prescriptions for the transition to nonperturbative trans¬ 
verse momenta 

We now focus on exploring the possibility of discriminating among different prescriptions 
that describe the functional dependence in bx of the nonperturbative Sudakov evolution 
factor (see Eqs. (3.20) and (3.21)) or the transition from the perturbative low— br domain 
to the nonperturbative high —br one (see Eqs. (3.18) and (3.19)). 

In Fig. 5, the normalized multiplicity of Eq. (5.2) is shown as a function of P 2 \ = 
zfq^ = (0.5) 2 g|i at the Belle scale Q 2 = 100 GeV 2 with the {bJp, 5 ^} prescription. Again, 
as in Fig. 3 there are two groups of uncertainty bands. The former one displays the results for 
the "fixed scale" evolution scheme in the standard notation, i.e. for { 6 max = 1.5, 92 = 0.18} 
(dot-dashed borders), {6 max = 1, 92 = 0.43} (dashed borders), and {6 max = 0.5, 52 = 0.68} 
(solid borders). The two additional light-gray bands correspond to the results with the 
" n\, scale" evolution scheme for {6 max = 1.5, g 2 = 0.18} (dot-dashed borders) and {6 max = 
0.5, g 2 = 0.68} (solid borders). So, also for the {b^ T , prescription we find the same 
ambiguity as for the {bj,, flnp} one i n Fig. 3: the overlap of the light-gray band with solid 
borders and of the band with dot-dashed borders indicates that two different evolution 
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Figure 5. The normalized multiplicity at Zi = Z2 = 0.5 as a function of PL = zWt = (0.5 )Wt 
at the Belle scale Q 2 = 100 GeV 2 for the "fixed scale" evolution scheme and with the {6^, 
prescription for the transition to the nonperturbative regime (see text). Notation and conventions 
for the uncertainty bands as in Fig. 3. 


schemes give similar results with different evolution parameters {6 m ax, <72}- Hence, we 
wonder if this similar trend suggests that it might not be possible to distinguish between 
the two schemes. Again, the possible way out is to look at the dependence of the results 
upon the fractional energy of the final hadrons. 




Figure 6. The normalized multiplicity at Z 2 = 0.5 as a function of P 2 ± = zfq f. at the Belle 
scale Q 2 = 100 GeV 2 for the ”fXb scale" evolution scheme and with the {btp, prescription for 
the transition to the nonperturbative regime (see text). Notation for the uncertainty bands as in 
previous figure. The additional light-gray bands with dot-dashed and solid borders are the result 
with the { 6 |,, < 7 }“} matching prescription for { 6 max = 1.5, 52 = 0.18} and {6 max = 0.5, <72 = 0.68}, 
respectively. Left panel for z\ = 0.5, right panel for Zi = 0.3. 


In Fig. 6 , the normalized multiplicity of Eq. (5.1) is shown as a function of P‘L = z\ 
at the Belle scale Q 2 = 100 GeV 2 for the " 7 ^ scale" evolution scheme. Also in this plot, 
there are two groups of uncertainty bands. A group displays the results for the { 6 ^, 
prescription in the standard notation, i.e. for {6 max = 1-5, <72 = 0.18} (dot-dashed borders), 
{5max = 1, #2 = 0.43} (dashed borders), and {6 m ax = 0.5, 52 = 0.68} (solid borders). The 
group of two light-gray bands correspond to the results with the {5^, t/J”)} prescription 
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for {& m ax = 1.5, <?2 = 0.18} (dot-dashed borders) and {6 max = 0.5, 52 = 0.68} (solid 
borders). If we focus on the left panel where calculations are performed at z\ = Z 2 = 0.5, 
the two bands with dot-dashed borders are substantially overlapped, thus reinforcing the 
suspect that it might not be possible to discriminate between the { bj,, <?}“} and {6^, ^njf} 
prescriptions. But if we now turn to the right panel, where the same calculation is performed 
at Z\ = 0.3, Z 2 = 0.5, we may hope to have a sufficiently small experimental error that 
discriminates between the two bands with dot-dashed borders. Unfortunately, the plot 
suggests also that this option seems possible only for the {6 max = 1.5, <72 = 0.18} case. And 
further explorations show that the same calculation, when performed in the "fixed scale" 
evolution scheme, produces more confused results. In summary, a combined study of the z 
and P^j_ dependencies in the normalized multiplicity might be able to discriminate among 
different prescriptions for the nonperturbative effects in the evolution only for a selected 
set of evolution parameters and schemes. 

5.4 Sensitivity to hadron fractional-energy dependence 

In the previous sections, we found that in several occasions only the combined study of the z 
and P 2 ± dependencies of the normalized multiplicity allows for discerning results obtained 
from different parametrizations and prescriptions in the description of nonperturbative ef¬ 
fects in the TMD evolution. This is not accidental. With the approximations adopted in 
this work, the main difference between the two considered evolution schemes lies in fact 
in the z dependence of the collinear fragmentation function d±, as it can be deduced by 
comparing Eqs. (3.27) and (3.28). 


M(z 1 ,z 2 ,pJ t ) / Mfz,,z,,0) MPi.Zz.pSt) / MfZj.Zz.O) 




Figure 7. The normalized multiplicity at Z 2 = 0.5 as a function of P 2 ± = zfq §, at the Belle 
scale Q 2 = 100 GeV 2 for the evolution parameters { 5 max = 1.5, 32 = 0.18} and with the {b^, g} 1 } 1 } 
prescription for the transition to the nonperturbative regime (see text). Uncertainty band with 
dot-dashed borders for z\ = 0.3, with dashed borders for Zj = 0.5, with solid borders for z\ = 0.7. 
The squared box with error bar corresponds to an experimental error of 7%. Left panel for the " 
scale" evolution scheme, right panel for the "fixed scale" one. 

The plots in Fig. 7 seem to confirm this finding. In the left panel, the normalized 
multiplicity of Eq. (5.1) is shown as a function of Pfj_ = zfq^ at the Belle scale Q 2 = 
100 GeV 2 for the "fib scale" evolution scheme, the {b^, gjjp} prescription, and the choice 
{6max = 1.5, (72 = 0.18}. The bands display results for the values z\ = 0.3, Z2 = 0.5 (band 
with dot-dashed borders), z\ = Z 2 = 0.5 (dashed borders), and z\ = 0.7, Z 2 = 0.5 (solid 
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borders). In the right panel, we show the results of the calculations performed in the same 
conditions but for the "fixed scale" evolution scheme. It is quite evident that the latter 
scheme produces P 2 j_ distributions that are systematically larger for any combination of 
{z\, £2}- This finding holds true also for other choices of the evolution parameters { 6 m ax> S 2 } 
and for the {btp, <?n°p } prescription. 

5.5 Sensitivity to the hard scale: from Belle to Bes-III 

All previous results have been obtained at the Belle scale of Q 2 = 100 GeV 2 . We may 
wonder what happens when reducing the "evolution path" to lower scales, like, e.g., the 
Bes-III scale Q 2 = 14.6 GeV 2 . 


M(z 1 ,z 2 ,Phr) / M(Z!,z 2 ,0) 



Figure 8. The normalized multiplicity of Eq. (5.2) as a function of P 2 ± = z 2 q^ = (0.5 ) 2 q^ at 
the Bes-III scale Q 2 = 14.6 GeV 2 for the "/ib scale" evolution scheme and with the {5^, <?)“} 
prescription for the transition to the nonperturbative regime (see text). Notation and conventions 
for the uncertainty bands as in Fig. 2. 

In Fig. 8, the normalized multiplicity of Eq. (5.2) is shown as a function of P 2 ^ = 
zfq^ = (0.5) 2 q|n in the same conditions and notation as in Fig. 2 but at the Bes-III 
scale Q 2 = 14.6 GeV 2 . By comparing these results with the ones in Fig. 2, we deduce 
that the net effect is a systematic enlargement of the uncertainty bands. This finding 
occurs also for other combinations of evolutions schemes and nonperturbative prescriptions. 
Hence, we deduce that working at the Bes-III scale is not useful if we want to discriminate 
among different evolution parameters {6 ma y, < 72 }, or between the {bj*, gjjp} and {6^, 5n°p S } 
prescriptions, or between the "fixed scale" and "scale" evolution schemes. 

However, we recall that each uncertainty band is the envelope of the 68% of 200 different 
curves, each one corresponding to a specific replica of the intrinsic parameters entering the 
Gaussian widths (Pf) a ^ h (z) of Eq. (3.1) for the bx distribution of the at the starting 
scale in the evolution. Then, we might envisage that the experimental error is sufficiently 
smaller than the band width such that it is able to discriminate some of the replicas, in 
order to narrow the uncertainty on the intrinsic parameters. In any case, this goal will be 
achieved only by performing additional more precise measurements of SIDIS multiplicities 
for different final hadron species and on different targets. 
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5.6 Sensitivity to partonic flavor 

The sensitivity to the nonperturbative intrinsic parameters, that describe the bx distribution 
of the TMD FF at the initial scale of evolution, is an important issue. The analysis of SIDIS 
multiplicities at low Q 2 suggests that some of these parameters are different for different 
flavors [27]. Hence, we expect that also the distribution in transverse momentum space of 
the evolved TMD FF will depend on the flavor of the fragmenting partons. However, the 
cross section in Eq. (2.9) mixes all flavors in the sum. Therefore, it is useful to define an 
observable that is well suited to explore the effect of flavor in the TMD evolution. 

In the following, we will show results for the P 2 ± distribution of ratios of normalized 
multiplicities corresponding to different final states: 


M hlh *(z u Z2, 4, y)/M hlh *(z u Z2, o, y) x 


M h M ( Zl ,z 2 , <?t, y)/M h ' lh ' 2 (*, z 2 ,0, y) 


-l 


(5.3) 


MOrV) / M(K + K~) 

1.6 

1.4 

1.2 

1.0 

MOfUr - ) / M(K + K") 

fixed scale - flavor indep. fb}, gj]p) 
b™,=1.5 
g 2 =0.18 

fixed scale - flavor dep. lb}. gUp) 

bmm=1.5 

g 2 =0.18 


0.8 

. --SJ.!! --- 

—- Zi=0.3, z 2 =0.5 


- Zi=0.3, z 2 =0.5 

— Zi=0.5, z 2 =0.5 

0.6 

- Zj=0.5, z 2 =0.5 

- Zi=0.7, z 2 =0.5 


- Zi=0.7, z 2 =0.5 


n a 



0.4 L- 1 —*—*—*—*— 1 —*—*—*—*— 1 —*—*—*—*— 1 —*—*—*—*—0.4 1 - 1 —*—*—*—*—*—*—*—*—*—*—*—*—*—*—*—*—*—*—*—<-* 
0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 


P^IGeV 2 ] P= T [CeV 2 l 


Figure 9. Tha ratio of normalized multiplicities in Eq. (5.3) between the {7r + 7r~} final state and 
the {I\ + K~} final state at z 2 = 0.5 and y = 0.2 as a function of P 2 ± = z 2 q^ at the Belle 
scale Q 2 = 100 GeV 2 for the "fixed scale" evolution scheme, for the evolution parameters {6 max = 
1.5, g 2 = 0.18}, and with the {K}, g}™} prescription for the transition to the nonperturbative regime 
(see text). Uncertainty bands with dot-dashed, dashed, and solid borders for z\ = 0.3, 0.5, 0.7, 
respectively. Left panel for flavor independent intrinsic parameters of input TMD FF, right panel 
for flavor dependent ones (see text). 

In Fig. 9, we show the ratio of Eq. (5.3) between the normalized multiplicity for {7r + 7r - } 
and the one for {K + K~} at z 2 = 0.5 and y = 0.2 as a function of P 2 ± = z\q^ at the Belle 
scale Q 2 = 100 GeV 2 for the "fixed scale" evolution scheme, for the evolution parameters 
{frmax = 1.5, g 2 = 0.18}, and with the {&!}, gjU} prescription for the transition to the 
nonperturbative regime. 

If we suppose to switch off the flavor dependence of the intrinsic parameters, the bx 
distribution of the TMD FF in Eq. (3.28) is controlled by the same Gaussian width (P'j) (z) 
for all channels. This feature remains valid when performing the Bessel transform to mo¬ 
mentum space, such that the distribution of the cross section can be factorized out of 
the flavor sum. Therefore, if we take the ratio of normalized multiplicities at the same z\ 
we expect the latter to be independent of P 2 ± = z 2 q^ ,. This is indeed the result displayed 
in the left panel of Fig. 9. It is a systematic feature of the "fixed scale" evolution scheme: 
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it holds true for other values of z\, as shown in the panel, but also for other combinations 
of nonperturbative evolution parameters and nonperturbative prescriptions. 

If we account for the flavor dependence of the Gaussian widths (Pj) q ~^ h (z), then the 
bx distribution is different for the {7r + 7r - } final state from the one for {K + I\~}, as it 
can be realized by inspecting Eqs. (4.29)-(4.35). Consequently, the ratio of normalized 
multiplicities has a specific P^ ± = z\qj* distribution that, of course, changes with z\. This 
is indeed the content of the right panel in Fig. 9: the uncertainty band of the 68% of 200 
replicas of Gaussian widths with dot-dashed borders corresponds to z\ = 0.3, the band with 
dashed borders to z\ = 0.5, the band with solid borders to Z\ = 0.7. 

Almost all the ratios are smaller than unity because in our approximations the frag¬ 
mentation into kaons has two favoured channels while the fragmentation into pions only 
one (see Eqs. (4.29) and (4.34)), and the P\ ± distribution of the fragmentation into kaons 
seems to be larger than the corresponding one for pions (see the analysis of Ref. [27]). In 
any case, we believe that the inspection of the Pf ± distribution of ratios of normalized 
multiplicities for different final hadrons produced in future e + e _ annihilation experiments 
is a useful tool to discriminate among different scenarios in TMD evolution. For example, if 
future data for this observable will lie well above unity, the "fixed scale" evolution scheme 
would be ruled out, independently of the flavor dependence of the intrinsic parameters in 
the TMD FF at the initial scale of evolution. 

In Fig. 10, in the two panels of the upper row we show the same ratio of normalized 
multiplicities in the same conditions and notation as in the previous figure but for the 
"fj,b scale" evolution scheme. The left panel still corresponds to the case when the flavor 
dependence of the intrinsic parameters is neglected. However, in the "n f, scale" scheme the 
bT distribution of the TMD FF is influenced also by the collinear part of the fragmentation 
function: the d\^ h in Eq. (3.27) is evaluated at the running scale which is related to bx 
via Eqs. (3.11), (3.18), (3.19). Hence, when performing the Bessel transform of D\ in the 
cross section, the resulting , distribution depends on the flavor of the fragmenting parton 
even if the intrinsic parameters do not. This "perturbative" flavor dependence, induced 
by RGE acting on the evolved collinear part of the TMD FF, mixes with the possible 
flavor dependence of the intrinsic parameters, making it rather difficult to disentangle the 
two effects. The left panel in the upper row shows the ratio of normalized multiplicities 
as a function of P\\ = z\q^ for three different values of z\. As in the previous figure, 
the band with dot-dashed borders corresponds to z\ = 0.3, the band with dashed borders 
to z\ = 0.5, and the band with solid borders to z\ = 0.7. Suprisingly, all the ratios are 
larger than unity. When including also the flavor dependence in the intrinsic parameters, 
the uncertainty bands become larger because there is a marked sensitivity to all possible 
replica values of the intrinsic parameters themselves. Again, as in the previous section we 
can argue that experimental data will have a sufficiently small error to discriminate among 
the various replicas. 

A further constraint can be achieved by considering a different combination of final state 
hadrons in the ratio of normalized multiplicities in Eq. (5.3). The lower panel in Fig. 10 
shows the results for the ratio between a {7r + 7r - } final state and a {ir + K~} final state 
when neglecting the flavor dependence of intrinsic parameters of the TMD FF at the initial 
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Figure 10. Upper panels: same as in previous figure but for the "^ scale" evolution scheme. 
Lower panel: the ratio between the normalized multiplicities M 71 77 ( 21,22 = 0.5 ,q^,y = 
0.2)/M v+77 ~ ( 21,22 = 0.5,0, y = 0.2) and M 77+K ~(z 1 ,z 2 = 0.5 ,q%,y = 0.2)/M n+K ~ (zi,z 2 = 
0.5,0, j/ = 0.2) as a function of P‘f ± = zfqj, at the Belle scale Q 2 = 100 GeV 2 in the same 
conditions and with the same notation as in the upper panels, but for flavor independent intrinsic 
parameters of input TMD FF (see text). 


scale. The notation and conventions are the same as in the other panels. All the ratios 
are now lower than unity. Hence, combining this result with the content of the upper left 
panel could represent a very selective test of the "^ scale" evolution scheme. In fact, when 
neglecting the flavor dependence of intrinsic parameters the Pf ± distribution of normalized 
multiplicities for the {7T + 7r - } final state should be larger than the one for {K + I\ ~} at any 
21 , while at the same time it should turn out narrower than the one for {tt + K~} at any 21 . 
Moreover, if future data for the {7r + 7r - } back-to-back production in e + e - annihilation will 
display a much narrower P ^ distribution than for the {K + K~} production, at least by 
20%, this will represent a further selective test for calculations performed in this evolution 
scheme, as it can be deduced by combining the results in the panels of the upper row. 

Finally, we notice that because of charge conjugation symmetry (see Sec. 4.1) we pre¬ 
dict that the ratio between normalized multiplicities leading to (7r + , K~) and ( 7 r~, K + ) 
final states should be equal to unity, irrespective of the choice of evolution schemes, non- 
perturbative evolution parameters and prescriptions. It would be interesting to cross-check 
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this prediction by measuring this ratio as a function of 

6 Conclusions 

In this paper, we consider the semi-inclusive production of two back-to-back hadrons in 
electron-positron annihilations. We study the transverse momentum distribution of such 
pairs of hadrons by observing the mismatch between their collinear momenta, and we focus 
on charge-separated combinations of pions and kaons. We conveniently define the multiplic¬ 
ities in electron-positron annihilations as the differential number of back-to-back pairs of 
hadrons produced per corresponding single-hadron production, in analogy to the definition 
of multiplicity in SIDIS process. In particular, we analyze the multiplicities normalized to 
the point of vanishing transverse momentum in order to extract clean and uncontaminated 
details on the transverse momentum dependence of the functions describing the fragmenta¬ 
tion process (transverse-momentum dependent fragmentation functions - TMD FFs). The 
normalized multiplicities are advantageous also because they turn out to be almost insen¬ 
sitive to the theoretical uncertainty related to the arbitrary choice of the renormalization 
scale. 

We consider electron-positron annihilations at large values of the center-of-mass (cm) 
energy, namely in the experimental conditions of the Belle and Bes-III experiments. We 
study how TMD FFs evolve with the hard scale. The input expression for TMD FFs is 
taken from a previous analysis of SIDIS multiplicities measured by Hermes at low energy, 
which is assumed as the starting scale. Since the hard scale in annihilation processes is 
much larger, we perform realistic tests on the sensitivity to various implementations of 
TMD evolution available in the literature. 

We find that within a specific evolution scheme the transverse momentum distribution 
of normalized multiplicities at the Belle scale can be very sensitive to the choice of the 
parameters describing the nonperturbative part of the evolution kernel. A hypothetical 
7% error in such data (compatible with the observed experimental error in collinear back- 
to-back emissions in electron-positron annihilations) could discriminate among different 
choices of parameters that are justified and adopted in the literature. 

But we observe also that at the same Belle scale different evolution schemes with 
different nonperturbative parameters can give overlapping transverse momentum distribu¬ 
tions. Our global results indicate that different evolution schemes can be discriminated only 
by considering the combined dependence of normalized multiplicities on both the transverse 
momentum and the fractional energy carried by the final hadrons. And this finding holds 
true (with some limitations) also for the purpose of discriminating among different pre¬ 
scriptions for describing the transition from nonperturbative to perturbative regimes in 
transverse momentum. 

The dependence on the fractional energy of the final hadrons is contained in the collinear 
part of the TMD FFs. Different evolution schemes produce different evolution effects also 
in the collinear fragmentation functions, which in turn emphasize the differences in the final 
transverse momentum distribution of evolved TMD FFs. The dependence on the fractional 
energy is contained also in the average squared transverse momenta that describe the width 
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of the input distribution of the TMD FFs at the starting scale. Therefore, by studying this 
dependence it may be possible to reduce the uncertainty on the intrinsic parameters that 
describe these input distributions. 

To this purpose, focusing on the normalized multiplicities at the Bes-III scale looks 
more promising. In fact, we observe that in stepping down from Belle to Bes-III scale the 
transverse momentum distributions of normalized multiplicities become much more sensi¬ 
tive to the details of the input distribution at the starting scale. The uncertainty in the 
determination of the intrinsic parameters needed to fit the Hermes SIDIS multiplicities 
reflects in a larger spread of normalized multiplicities as functions of transverse momentum. 
At the Bes-III scale, a hypothetical experimental error of 7% does not discriminate among 
results coming from different nonperturbative evolution parameters or from different evo¬ 
lution schemes. But within a specific choice of evolution scheme it can discriminate among 
results that come from different values of the intrinsic parameters. 

The Hermes results al low energy show significant differences between SIDIS mul¬ 
tiplicities for final-state pions and kaons. Hence, these data were fitted using transverse 
momentum distributions for the input TMD FFs that contain flavor dependent parameters. 
Here, we explore also how the final results for normalized multiplicities at Belle and Bes- 
III scales are sensitive to the details of this flavor dependence at the starting scale. In doing 
so, we find that the most convenient observable is represented by the ratio of normalized 
multiplicities for different final hadron species, particularly at the Belle scale. 

The most striking evidence is for evolution schemes where the flavor dependence is 
strictly localized only in the intrinsic parameters of the input TMD FFs at the starting scale. 
If we switch off such flavor dependence, the transverse momentum distribution of normalized 
multiplicities is always the same, irrespective of the species of final hadrons. So, if we select 
for example pions and kaons, the ratio of the corresponding normalized multiplicities is 
constant and equal to unity. If the flavor dependence of the intrinsic parameters is switched 
on, then the ratio deviates to values (mostly) lower than unity, in agreement with general 
expectations that kaons have a larger distribution in transverse momentum. 

The situation is more confused for evolution schemes where the flavor dependence 
is indirectly contained also in the initial conditions of the evolution equations through 
the (flavor dependent) collinear part of the fragmentation functions. In this case, this 
effect mixes up with the flavor dependence contained in the intrinsic transverse momentum 
distribution, and it is difficult to disentangle one from the other. At variance with the 
previous class of evolution schemes, in this case the ratio of normalized multiplicities for 
pions with respect to kaons turns out to be (mostly) larger than unity. Fortunately, more 
selective criteria are offered by considering a variety of species of final hadrons. If we 
consider ratios of normalized multiplicities for pions with respect to mixed pion-kaon pairs, 
the results are (mostly) lower than unity. By combining the results for various final states 
all together, one would hope to constrain the arbitrary ingredients of TMD FFs as much as 
possible. 

We conclude by stressing that all the results and remarks above refer to the unpolarized 
TMD FFs that describe the fragmentation of an unpolarized parton into an unpolarized 
hadron. However, this function is an essential ingredient in all the (spin) azimuthal asym- 
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metries extracted in hard processes like electron-positron annihilation, hadronic collision, 
and SIDIS. Hence, a better control on the transverse momentum dependence of unpolarized 
TMD FFs implies also a better knowledge of polarized TMD FFs as well as of (un)polarized 
TMD parton distributions. For this reason, we are looking forward to a multidimensional 
analysis of data accumulated by the Belle and Bes-III collaborations, possibly including 
a study of normalized multiplicities for various hadron species as suggested in this work. 
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